Bright Solitary-Matter- Wave Collisions in a Harmonic Trap: Regimes of Soliton-like Behaviour 
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Systems of solitary-waves in the ID Gross-Pitaevskii equation, which models a trapped atomic Bose-Einstein 
condensate, are investigated theoretically. To analyse the soliton-like nature of these solitary- waves, a particle 
analogy for the solitary- waves is formulated. Exact soliton solutions exist in the absence of an external trapping 
potential, which behave in a particle-like manner, and we find the particle analogy we employ to be a good 
model also when a harmonic trapping potential is present. In the case of two solitons, the particle model is 
integrable, and the dynamics are completely regular. The extension to three particles supports chaotic regimes. 
The agreement between the particle model and the wave dynamics remains good even in chaotic regimes. In the 
case of a system of two solitary waves of equal norm, the solitons are shown to retain their phase diff'erence for 
repeated collisions. This implies that soliton-like regimes may be found in 3D geometries where solitary waves 
can be made to repeatedly collide out of phase, stabilising the condensate against collapse. 

PACS numbers: 03.75.Lm, 05.45.-a, 45.50.Tn 



I. INTRODUCTION 

Solitary- waves may be found in solutions to nonlinear wave 
equations where the nonlinearity counteracts the dispersion of 
a wave-packet such that it retains its form as it propagates. 
Solitons are solitary-waves that emerge unscathed from colli- 
sions with each other, up to shifts in position and phase; this 
behaviour is reminiscent of particle behaviour, motivating the 
particle-like name soliton. This distinction is an important 
one, although in practice the names soliton and solitary-wave 
are commonly interchanged. "Classic" solitons, in this sense, 
are to be found in integrable nonlinear wave-equations, such 
as the Korteweg-de Vries equation, the sine-Gordon equation, 
and the one-dimensional nonlinear Schrodinger equation. The 
solitons' ability to re-emerge after collisions is due to the fact 
that their dynamics are strongly constrained by conservation 
laws associated with the wave-equations' integrability 1 1]. 

Solitons and solitary-waves are topics of keen interest in 
the atomic Bose-Einstein condensate (BEC) community. This 
is because low-temperature BEC dynamics are frequently de- 
scribed to a good approximation by the Gross-Pitaevskii equa- 
tion (GPE) H i, 11, a 3D nonlinear wave equation. For 
regimes where the atoms are confined in the radial direction 
by a tight trapping potential, the 3D GPE reduces approx- 
imately to a ID equation (the so-called ID GPE). The ho- 
mogeneous ID GPE is simply the ID nonlinear Schrodinger 
equation, which can be solved by the inverse scattering trans- 
form, and yields bright soliton solutions when the nonlinearity 
is attractive | LlS]. At sufi&ciently low temperatures the inter- 
atomic scattering can be largely described through a single 
parameter, the 5-- wave scattering length. In this context, an 
attractive nonlinearity arises from a negative ^-wave scatter- 
ing length, which may be intrinsic, or which may be induced 
by exploiting a Feshbach resonance to tune the inter-atomic 
interactions |6, 7]. As well as describing BEC under tight 
transverse confinement, the ID nonlinear Schrodinger equa- 
tion is also used to describe nonlinear optical systems @] • 
These systems provide a useful analogue of BEC under tight 
transverse confinement, and we will frequently refer to work 
on nonlinear optics in this paper. 



Experiments involving BECs composed of attractively in- 
teracting atoms have been performed in ID geometries, re- 
sulting in the observation of single |10] and multiple bright 
solitary- waves iUl [13, IHl]. In the experiments with multi- 
ple solitary- waves, the BEC was trapped in the axial direction 
by a (relatively weak) harmonic confining potential in addi- 
tion to the radial confinement. The addition of an axial har- 
monic potential acts to break the integrability of the ID GPE, 
meaning that we no longer have exact soliton solutions. In the 
experiment by Strecker et al. classic soliton-like be- 

haviour (where the solitary-waves collide and reform up to 
shifts in phase and position) was not observed, but rather, 
trains of solitary-waves which are continuously repelled by 
each other. The dynamics of solitary-wave trains both in BEC 
and nonlinear optics have been the topic of extensive mod- 
elingusing a variational method lll4|] . numerical simulations 
IBlIm [O, fls", ^20], a Toda lattice approach, tlU, a par- 
ticle model 1 22] (quite distinct to that presented in this pa- 
per), analysis using the inverse- scattering transform ['23'] and 
by using a perturbation approach JH, |25l ^6] . These treat- 
ments model regimes where the solitary-waves are never well 
separated, where it has been found that the solitary-waves do 
not collide with each other and re-form, but interact with each 
other by attractive and repulsive forces, depending on their 
relative phase. Motivated by the observation of such soliton 
trains, a "soliton laser" has been proposed ll27h . A review ar- 
ticle on BEC solitons addresses some of this work in more 
detail ||28D. 

As opposed to solitary-wave trains, we investigate whether 
classic soliton-like behaviour, i.e., colliding and reforming of 
distinct, localized wave packets up to shifts in phase and posi- 
tion, is possible in the ID GPE with a harmonic potential. In a 
previous work ll29ll we found regimes where such behaviour is 
quite pronounced. This behaviour was also seen in work done 
in similar nonlinear optical settings [HO, [HI [H, [lH]. In this 
paper we further our investigation into soliton-like behaviour; 
in particular we explore the bounds within which the solitary- 
waves can still be expected to behave as solitons. To this end, 
we use a particle model introduced in our previous work |29|], 
adapted from a model developed for use in nonlinear optics 
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[31I ^5] . We show that soHton-Hke behaviour is possi- 
ble in the ID GPE with a harmonic potential, provided that 
the solitary-waves collide with large enough relative velocity 
such that the collisions occur during a short timescale com- 
pared with the period of the axial trapping potential. This typ e 
of behaviour has recently been experimentally observed , 
and provides an exciting prospect for future experiments to 
probe the dynamics in more detail. 

In the case of three solitons, we find regimes of regular and 
chaotic dynamics. In particular, chaotic solutions to the GPE 
are expected to coincide with more rapid condensate depletion 
than in otherwise similar regular solutions [34]; indeed this 
has been seen in theoretical studies of several systems [i35[ [36[ 
[STllIslllil]. This provides an additional motivation to identify 
regimes of regular and chaotic soliton dynamics in the GPE. 

In more realistic models for BECs, the integrability of the 
nonlinear wave equation is also broken by residual 3D ef- 
fects. These eff'ects cause the soliton collisions to be inelastic; 
specifically, there is particle exchange between the solitons 
accompanied by changes in their outgoing velocities ||4Q[]. A 
reduction from 3D to non-integrable ID equations [4r, ^2], 
more sophisticated than the ID GPE, confirms this result [43]. 
This type of behaviour is common in other non-integrable 
Schrodinger- type equations: 10, SI, H?!]. Bose-Einstein 
condensates with attractive interatomic interactions are prone 
to collapse if the particle density becomes too high [11]. Fully 
3D GPE simulations show that in-phase collisions between 
solitons, during which the particle density becomes large, can 
cause collapse of the condensate, [111 SO, • In this paper, 
as well as identifying regimes in which soliton-like behaviour 
can still occur despite the mild breaking of integrability by 
the harmonic potential, we also briefly discuss regimes where 
solitons are expected to survive 3D integrability breaking. 

The layout of the paper is as follows: in section II we in- 
troduce the model equations and reiterate the soliton solution 
to the homogeneous ID GPE; in section III we identify one, 
two and three solitary- wave solutions to the ID GPE with a 
harmonic potential and introduce a particle model to test their 
soliton nature; in section IV we present our conclusions. 



g3D = Anfp-alm and m is the particle mass of the species. De- 
pending on the species, a may be positive or negative, corre- 
sponding to an efl'ective repulsive or attractive interaction. By 
exploiting a Feshbach resonance, it may also be tuned using 
an external magnetic field 11,01]. 

The many-body Hamiltonian for the system in second- 
quantized form is then given by: 



^(r). 



where 



^ = -^ + Vext(r) 
2m 



(1) 



(2) 



is the Hamiltonian for a single particle in the external trap- 
ping potential. If Bose-Einstein condensation has occurred, 
we may define the condensate mode as the eigenfunction ^(r) 
of the single-body density matrix p(r, r') = (^^(r')^(r)) with 
the largest eigenvalue 15^ in a. We are then free to par- 
tition the fieldoperator into condensate and non-condensate 
parts 



(3) 



where a annihilates a particle in mode ^(r), and ^^'(r) is the 
field operator for modes orthogonal to the condensate. 



2. 3D classical field 

In the case of a trapped, almost fully Bose-condensed dilute 
atomic gas, the dynamics of the condensate mode ^(r) are 
largely governed by the following GPE 1531 154|]: 



ot 



2m 



+ Vext(r) + g3DA^|^(r,0r 



^(r), (4) 



II. BACKGROUND 
A. Model system 

1. Quantum field 



where N is the total number of particles in the condensate and 
^(r) is normalised to one. 

We now consider a cylindrically symmetric (cigar- shaped) 
harmonic trapping potential: 



(5) 



In the case of an atomic Bose gas, the dynamics of the sys- 
tem may be described, in the Heisenberg picture, by the time- 
evolution of the bosonic field operator ^(r). At the tempera- 
tures typically encountered in atomic BEG experiments (nK), 
for sufl&ciently dilute Bose gases, atomic interactions are dom- 
inated by low energy two-body collisions. In this case, the 
atom-atom interactions are characterised by one parameter: 
the 5"- wave scattering length, a. Moreover, we may generally 
replace the true interaction potential by an eff'ective contact in- 
teraction, subject to an appropriate renormalization procedure 
EIEIISIJ. Hence, we let yint(r,rO = g^j^dir - rO where 



where cOx ^ (j^^r- We also explicitly assume a < (attractive 
inter-particle interactions), and determine that Eq. © reduces 
to the following ID equation (see appendix (A) for a deriva- 
tion): 

ij^H^) = - 2 + -^H^) - I^AWI VW, (6) 

where x is measured in units of ff-fmlgwlN and t in units of 
/mlgiD^-N^, with giD = 2fL(ji)rCi and co equal to the axial 
frequency co^ in our units of inverse time im\gix)^N^ jli?). 
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3. Linear instability and chaos 



It can be shown that Hnear instabiUties in the GPE di- 
rectly imply, in both Bogoliubov 1 56] and equivalent number- 
conserving linearised approaches |53], that the population of 
the non-condensate component may rapidly become signifi- 
cant. For this reason we expect regimes where the GPE dy- 
namics are chaotic to coincide with rapid depletion of the con- 
densate. The GPE is a norm-conserving equation, so the de- 
pletion will not show up in the GPE dynamics; however, we 
may use chaos in GPE dynamics as an indicator of depletion 
of the condensate mode in a realistic system. This motivates 
the identification of chaotic trajectories in the GPE dynamics, 
as we discussed in the introduction (SeclJ). 



B. Soliton solution to homogeneous GPE 

In the case of no axial potential [equivalent to setting co = {) 
in Eq. Q], it is possible to find exact solutions of Eq. Q 
111]. A straightforward interpretation is in terms of a scattering 
problem [1]; in the limit ^ ^ -oo, the solutions take the form 
of an arbitrary number of well separated (incoming) solitons: 



O/jc, = 27/ysech ^rij{x - qj)\ e 



ivj(x-qj) i(2T]j +vj/2)t iaoj 



(7) 



Here qj = Vjt-\-xoj is the position of the peak of the 7th soliton; 
xoj is the peak position at ^ = 0; aoj - vjXqj is the phase for 
a single soliton (i.e., in the absence of collisions with other 
solitons) at X = 0, ^ = 0; Vy is the soliton velocity and T]j 
gives the relative size of the soliton. Our normalisation condi- 
tion implies ^Vj = 1 ^ where A^^ is the number of solitons 
present. 

The solitons come together and collide, during which time 
the form of the solution is complicated and solitons are not 
individually defined. However, as ^ cx), the outgoing soli- 
tons re-emerge from the collisions unscathed, taking the same 
asymptotic form [Eq. d?])], up to shifts in position and phase: 

^ + ^^j ^nd aoj aoj -h Scpj, where the position shift 
Sxj and phase-shift Scpj of the 7th soliton are given by IlllSt]: 



2r]Sxj -h iScpj = +2 In 



Vk + i2(r]j -\- rjk) 



Vj-Vk-\-i2(r]j-r]k) 



(8) 



The positive sign applies if the soliton is on the left prior to 
the collision with the ^th soliton (vj > Vk), otherwise the neg- 
ative sign applies. Note that these shifts are dependent on the 
solitons' initial speeds vj, and eff'ective masses T]j only, not on 
their relative phase. 



III. SOLITON DYNAMICS WITH A HARMONIC 
EXTERNAL POTENTIAL 

A. Single soliton 

As shown in appendix B, for any solution to the ID har- 
monic GPE, there exist other solutions with the same density 



0.25 



0.2 



1^0.15 



0.1 
0.05 



-10 




10 



20 



FIG. 1: (Colour online). Ground state solution for harmonically 
trapped soliton of 5000 particles (solid line). Corresponding soli- 
ton solution to the homogeneous equation (dot- dashed line), which 
is used as an ansatz in the particle model. The ground state of the lin- 
ear Schrodinger equation (dashed line) is given for comparison. The 
parameters of the system are taken to be similar to those of a recent 
experiment [13] the axial trapping frequency is 10/2:7r Hz, the radial 
trap frequency is 800/2:/r Hz, atomic mass and scattering length of 
^Li. A unit of x is hence equal to 7.19 X 10~^ m. 



profile, undergoing arbitrary amplitude harmonic oscillations 
at the trap frequency. In particular, for any stationary solution, 
there exist corresponding solutions with the same density pro- 
file, which oscillate with the trap frequency but remain other- 
wise unchanged. Hence, a single bright soliton in a harmonic 
trap experiences an overall simple harmonic motion without 
any manifestation of internal dynamics in the soliton 's den- 
sity profile. 

The density profile and phase behaviour of a single soli- 
ton can be found by first considering the form of a stationary 
soliton, and then inferring the behaviour of the oscillating ver- 
sion. The stationary soliton will be a solution to the eigenvalue 
problem: 



1 ^2 



2 2 



i/f(x)-mx)U(x)=/n/i(x\ (9) 



and will have the form: i/^(x) = u(x)exp{-i[jit -\- S(0)]}, 
where u(x) is a real valued function and the real- valued num- 
ber S (0) is an initial phase. We expect the single stationary 
soliton solution to be the metastable "ground" state of the sys- 
tem which may be determined numerically, for 
example by propagating Eq. © in imaginary time [60, 61]. 
The numerically determined density u(x)^ for a parameter 
regime consistent with the "^Li experiments of Strecker et al. 
is shown in figure [T] and is compared to a bright soliton so- 
lution of the homogeneous ID GPE, and to the ground state 
of the ID linear Schrodinger equation with harmonic poten- 
tial. As expected, the solution of the ID GPE with a har- 
monic potential is spatially slightly compressed compared to 
the bright soliton solution of the homogeneous GPE. These 
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two solutions, however, are quite similar (and can be made 
more similar as cd is progressively reduced), and are quite 
distinct from the Gaussian solution produced by the linear 
Schrodinger equation. We will exploit this similarity later in 
the paper. 

As shown by the treatment in appendix |Bl an oscillating 
soli ton solution takes the form: 

il/{x, t) = u{x- <x(0)] exp (/ [-jdt + {p{t))x - S (t)]) , (10) 

where {x(t)) = xq cos(ojt) + (po/oj) sin(6i;0 is the position ex- 
pectation value of the atomic ensemble, {p(t)) = po cos(ojt) - 
ojxq sin(6i;0 is the momentum expectation value, i.e.. 



-f 



{x) = I dxi//*(x)xi//(x) 



and 



(p) = -i I dxi^*(x) 



^dil/{x) 
dx 



(11) 



(12) 



po and xq are the initial position and momentum expectation 
values, respectively, and 

c/A { 2 sin(2a)f) xqPq ^oPo^„,„, 

Sit) = [P,-—] + — cos(2..0 - — + 5(0). 

(13) 

Removing the nonlinearity reproduces the result for a coher- 
ent state. When the nonlinearity is present, however, the sta- 
tionary eigenvalue, yu, is dependent on the norm of the soliton, 
unlike the case of the linear Schrodinger equation. 



B. Two solitons 

1. Overview 

In this section, we present a simple model of multiple 
trapped solitons, treating each of the solitons as a classical 
particle. We explore the case of two harmonically trapped 
solitons, and present results comparing the trajectories in the 
particle model with simulations of the wave dynamics in the 
GPE. These results allow us to determine the range of initial 
conditions for which the particle model is a good description 
of the system. 



2. Particle model 

Recall from Sec.|ll]that, in a homogeneous system, the tra- 
jectories of solitons emerging from collisions with each other 
are independent of the relative phase of the incoming solitons. 
The only effect of the relative phase of the solitons is on the 
form of the wavefunction (peak or trough) during the colli- 
sion. This is illustrated in figure [S] The phase-independence 
of the solitons' incoming and outgoing trajectories allows a 
model to be formulated that treats the solitons as classical par- 
ticles, each with only the positional degree of freedom (rather 



-(qj-qk) 



(14) 



than position and phase degrees of freedom used, for example, 
in f2^). This model was introduced by Scharf and Bishop in 
the context of nonlinear optics (13, HH [12], which we have 
adapted for the purpose of modeling a quasi- ID harmonically 
trapped BEC [29]. 

To construct the particle model, we first consider the ho- 
mogeneous solution before introducing the eff'ects of the har- 
monic trap. Following the approach in [ 62], one can derive an 
eff'ective inter-soliton potential (see appendix iDl): 

^ r 2rijrik 

y(qj - qk) = -2r]jr]k(r]j + 7]k)sQch 

which treats the solitons as particles of position qj and eff'ec- 
tive mass T]j, the parameters used to describe the bright soliton 
solutions of Eq. (|7]). This potential reproduces the asymp- 
totic position shifts [Eq. ([8])] in the homogeneous GPE for the 
outgoing particle trajectories, i.e., the position shifts as the 
solitons become infinitely far apart. It yields accurate results 
when 2\t]i - 7/2I «c |vi - V2I, where vi and V2 are the soliton ve- 
locities, and which gives a lower limit for the relative velocity 
for which the particle model is applicable. 

Figure [2l shows the particle trajectories predicted by our 
model interaction potential [Eq.[T4l superimposed on the den- 
sity profile dynamics predicted by solution of the homoge- 
neous ID GPE. When modelling BEC dynamics, an upper 
limit to the solitons' relative velocity is also imposed, because 
the contact-interaction potential between atoms, used to de- 
rive the GPE, assumes low energy inter-atomic collisions, and 
may not be ap plic able to condensates with high relative ap- 
proach speeds ll63|] . Fortunately, recent experiments (0,13 
show that solitons are generated with similar sizes, such that 
their velocities may easily fall within our model's range of 
validity. 

In Sec, nil Al we showed that independent trapped solitons 
oscillate harmonically. A more general method for deriving 
the approximate motion of solitons in an external pot ential 
of arbitrary form was found by Scharf and Bishop |32], and 
is outlined in appendix [C] for completeness. To combine the 
eff'ects of the external potential and the soliton collisions, we 
use the homogeneous solution [Eq. (|7])] as an ansatz, so that 
the solitons are still characterised by the parameters qj and rjj. 
Figure [T] shows this to be a reasonable approximation. The 
following Hamiltonian: 



H 



(15) 



reproduces the harmonic motion of the solitons, keeping the 
interpretation of r]j as eff'ective masses (see appendix [C])- We 
assume that the soliton- soliton interactions are not aff'ected by 
the introduction of the (relatively loose) harmonic trap and 
construct the full Hamiltonian by adding in the contributions 
from the interaction potentials: 



H 



^ 27] j7]k(7] j + 7]k)SQCh^ 



l<j<k<Ns 



(qj-qk) 



(16) 
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FIG. 2: (colour online). Two soliton collision taking place when (a) A0coi = 0, (b) A0coi = 7t/2, (c) A0coi = and (d) A0coi = 3;r/2. The 
parameters of the system are taken to be similar to those of a recent experiment 1 13] the axial trapping frequency is 10/2:7r Hz, the radial trap 
frequency is 800/2;7r Hz, atomic mass and scattering length of ^Li, and 5000 particles per soliton). The unit of x is then equal to 3.6 fim, and a 
unit of ^ to 1.4 ms. 



where A^^ is the number of solitons. This approach is expected 
to be valid for regimes when the timescale of the soliton- 
soliton collisions is much less than the period of the harmonic 
trap, such that the effects of the harmonic trap are negligible 
during the collisions. The limits of this approach are further 
explored in Secs. lIIIB 3l and lIIIB 41 

In the case of two solitons (A^^ = 2), it is useful to define 
the following independent coordinates: the centre-of-mass po- 
sition Q '= (jiiqi + ^2^2)/(^i + Vi) ^iid the relative position 
q := qi - q2- The Hamiltonian [Eq. ([T6l) 1 then takes the form: 

ti=——- + —{m+m)Q^ 
2(m + 7]2) 2 



7/1+7/2 2 ^' 



27/17/2 



mm ^2 



2 7/1+7/2 



(17) 



- 27/17/2(7/1 + 7/2)sech 



27/17/2 

m +m 



where P = /7i + /?2 is the momentum canonically conjugate to 
2, and p = (r]2Pi - mP2)Km + m) the momentum conjugate 
to q. The Hamiltonian is now clearly seperable into two parts: 
the centre-of-mass energy E (dependent on P and Q only), and 
the interaction energy e (dependent on p and q only). There 
are thus two independent constants of the motion, E and 6, 
as many as there are degrees of freedom. Hence, the particle 



model for two solitons is integrable and the dynamics must 
be completely regular |65l] . In the case where the solitons 
have identical eff'ective masses (7/1 = 7/2 := 7/), the Hamilto- 
nian [Eq. ([T7l) 1 reduces to 



2 22 
2^2 . . T 



47/ 7/ 



-47/^sech2(7/^). (18) 



Figure [3] shows four Poincare surfaces of section (or 
Poincare sections) for the two particle system; sections of pi 
versus qi are shown for diff'erent surfaces of diff'erent total en- 
ergy 164, 65]. These Poincare sections demonstrate the regular 
behaviour of the integrable two particle system, as all trajec- 
tories lie on invariant tori in the phase space of the system. 
There are two distinct regimes observable in these Poincare 
sections. In the lower regions of the sections, the centre-of- 
mass energy, E, is large and positive; in this case the inter- 
action energy, 6, has a large negative contribution from the 
interaction term, and the solitons interact strongly. It is seen 
in section HTlB 3l that in this regime, there is rapid energy ex- 
change between the solitons, such that the soliton with lower 
amplitude oscillations is driven by the other soliton, which it- 
self becomes damped. In the upper regions of the sections, E 
is less positive, and hence e is less negative, so the energy ex- 
change between the solitons occurs over a longer time period. 
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FIG. 3: (Colour online) Poincare sections for the two-soliton system corresponding to the momentum pi and position qi of one soliton, while 
the other soliton has coordinates q2 = 0, /?2 < 0. The value of centre-of-mass energy, E, in the is given by the colour scale, (a) Total energy 
H = 5 X 10""^; (b) // ^ 5.6 X 10~^, the star corresponds to the trajectory in figure IH (c) // ^ 8.1 x 10~^, the upper trajectory correspond to that 
in figure|5] the lower to that in figure[6j (d) H ^ 2.2x 10~^, the star corresponds to the trajectory in figure|7] The figures correspond to regimes 
where the solitons have equal eff'ective masses, the axial trapping frequency is 10 /In Hz, and the other parameters (radial trap frequency of 
800/2:7r Hz, atomic species mass and scattering length of ^Li, and 5000 particles per soliton) are comparable to those in recent experiment iflSH . 



3. Presentation of results comparing GPE to particle evolutions 

In the Poincare sections of Fig. [3l we highlighted a num- 
ber of trajectories in white. These trajectories are plotted in 
position space as a function of time, overlaying density plots 
of corresponding ID GPE solutions, in Figs.|4]-[7l We do this 
to test how accurately the particle model represents the GPE 
dynamics for a range of initial conditions. 

In particular, the wave-dynamics involve a phase variable 
not accounted for in the particle model. As shown in Fig. [2l 
the phase diflTerence between the solitons has an observable 
effect on the solitons' form during collisons, although as the 
solitons tend asymptotically apart, the solitons' density dy- 
namics are insensitive to this. Two solitons with equal norms 
in a harmonic trap have the same collisional form for all sub- 
sequent collisions; i.e., two solitons initially colliding with a 
phase difference 0coib will have this phase difference for all 
subsequent collisions (see appendix [E]). This property allows 



results of GPE simulations with repeated in-phase and n out- 
of-phase collisions to be compared for any trajectory in the 
two-particle model. Regimes with phase differences between 
zero and tt are not considered here, but will generally be ex- 
pected to display behaviour intermediate between that of the 
zero and tt cases. 

The trajectories displayed in Fig. |4] correspond to the fixed 
point (marked by a white star) in the Poincare section of Fig. 
^h). We observe significantly better agreement between the 
particle model and the wave dynamics for the in-phase case 
[Fig-ffla)] than for the out-of-phase case [Fig.Ub)], where an 
obvious discrepancy between the particle and wave dynamics 
gradually accumulates. Subsequent collisions consistently oc- 
cur slightly earlier in the particle model [dynamics induced by 
the Hamiltonian ([TSt i than for the solitons propagated by the 
GPE [Eq. ©], so that by t = 5000 quite a noticeable shift of 
the particle dynamics has taken place. For the zero phase case, 
a similar systematic discrepancy is just observable at large t, in 
Fig.lHa); however, this discrepancy is very small and not read- 
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FIG. 4: (colour online). Trajectories in the particle model (lines) plotted over density distributions predicted by ID GPE dynamics, corre- 
sponding to the trajectory marked on figure Hb). The relative phase of the solitons in the wave dynamics is zero in figure (a), and n in figure 
(b). The figures correspond to regimes where the solitons have equal eff'ective masses, the axial trapping frequency is 10/2;7r Hz, and the 
other parameters (radial trap frequency of 800/2;7r Hz, atomic species mass and scattering length of ^Li, and 5000 particles per soliton) are 
comparable to those in recent experiment 1 13]. The unit of x is then equal to 3.6 jim, and a unit of t to 1.4 ms. 



ily apparent for most of the evolution. We note that t = 5000 
corresponds to 7 seconds using the parameters of Ref. (i^ , 
i.e., substantially longer than most typical experimental time 
scales. 

Figure[5] shows trajectories equivalent to those of Fig.|4l but 
with an additional centre-of-mass motion. This corresponds 
to the upper of the two trajectories marked in white in the 
Poincare section shown in Fig.[3tc), where there is gradual en- 
ergy exchange between the solitons. We see the same discrep- 
ancy between particle and wave dynamics as was observed in 
Fig- in This is to be expected, because, as shown in appendix 
IHl incorporating a centre-of-mass oscillation into a solution 
of the harmonic GPE simply causes the wave-function's den- 
sity profile to oscillate, without any further effect on its overall 
evolution. 

The trajectories in Fig. [6l correspond to the lower trajec- 
tory illustrated in the Poincare section of Fig. [Sjc), where 
there is rapid energy exchange between the solitons. In this 
regime, good agreement is found in the case of in-phase col- 
lisions. For the case of out-of-phase collisions, there is some 
degree of agreement, but there are obvious discrepancies. No- 



tably, the solitons in the GPE appear to remain a fixed dis- 
tance apart, whereas the particle trajectories necessarily show 
repeated collisions. 

The dynamics shown in Fig. [7] correspond to the fixed point 
(marked by a white star) in the Poincare section of Fig.[3td). 
In this regime, the solitons collide with a significantly faster 
relative speed than in Figs. |4]-[6l and the collision time is ac- 
cordingly shorter compared to the trap period. In this regime 
it is clear that the agreement between the particle model and 
the wave dynamics is much better than in the previous cases 
(figures |4] to O, with the relative phase having much less of an 
observable influence. 



4. Discussion of results comparing GPE to particle evolutions 

To explain the observed discrepancies between correspond- 
ing particle and wave evolutions, particularly when the soli- 
tons collide out-of-phase (see Figs. |4land[5]), we must consider 
the assumptions made while composing the eflTective particle 
Hamiltonian [Eq. ([T6b l . When constructing this Hamiltonian, 
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FIG. 5: (colour online). Trajectories in the particle model (lines) plotted over density distributions predicted by ID GPE dynamics. The 
trajectories correspond to those given in figured but with additional centre-of-mass displacements. The trajectories also correspond to the 
upper trajectory marked on figure [3c). The relative phase of the solitons in the wave dynamics is zero in figure (a), and n in figure (b). The 
solitons have equal eff'ective masses, the axial trapping frequency is 10/2;r Hz, and the other parameters (radial trap frequency of 800/2;7r Hz, 
atomic species mass and scattering length of ^Li, and 5000 particles per soliton) are comparable to those in recent experiment 1 13]. The unit 
of X is then equal to 3.6 yum, and a unit of tio 1.4 ms. 



we stated that the collision time should be small compared 
with the trap period. To characterise the timescale of a colli- 
sion, it is instructive to consider how rapidly the trajectories 
in the particle model converge, to those deduced from motion 
of the soliton peaks in the GPE, after a collision. 

We consider again the evolutions depicted in Fig. |2ta) and 
[3b). Figure [8] shows the differences in position between the 
peak of the left-hand soliton evolved by the GPE, and the cor- 
responding trajectory in the particle model during a collision. 
As there is no external potential, the particle trajectories are 
asymptotically exact as ^ ^ ±oo. It is clear that, in the case of 
an in-phase collision, the convergence of the particle and wave 
trajectories is more rapid than for a tt out-of-phase collision. 

In a harmonic trap, subsequent to a collision, two solitons 
can only move a finite distance apart (i.e., not asymptotically 
far) before moving together and colliding once more. The 
dynamics of solitons colliding with a n phase diflference are 
therefore not expected to agree as well with the effective par- 
ticle model dynamics, compared with solitons colliding in- 
phase. This eflTect builds up over time because the phase dif- 



ference is preserved for repeated collisions. Figure [8] shows 
that, in the n out-of-phase case, compared to the particle tra- 
jectory, the trajectory of the soliton peak tends to be further 
away from the point of collision, with the two trajectories con- 
verging asymptotically as t tends to a time infinitely before 
and after the time of collision. Due to the harmonic poten- 
tial, the particle trajectories consequently start their return to 
the centre of the trap at points closer to the centre, compared 
to the 71 out-of-phase solitons. This explains why subsequent 
collisions take place earlier for the particle trajectories than 
predicted by GPE dynamics, as observed in Fig. [H For large 
approach speeds (Fig. [3), however, the collision time is suf- 
ficiently small for both in-phase and out-of-phase cases, such 
that any discrepancy between the predictions of the GPE and 
the efifective particle model is too small to be observed. 

Figure [6] represents a regime where the solitons do not sep- 
arate well between collisions. The particle model might not 
be expected to apply well to such regimes; nevertheless, rea- 
sonably good agreement is observed in Fig. [3a), with some- 
what less good agreement (as expected) for the tt out-of-phase 
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FIG. 6: (colour online). Trajectories in the particle model (lines) plotted over density distributions predicted by ID GPE dynamics, corre- 
sponding to the trajectory marked on figure Hd). The relative phase of the solitons in the wave dynamics is zero in figure (a), and tt in figure 
(b). The solitons have equal eff'ective masses, the axial trapping frequency is 10/2;r Hz, and the other parameters (radial trap frequency of 
800/2;7r Hz, atomic species mass and scattering length of ^Li, and 5000 particles per soliton) are comparable to those in recent experiment [ 13ll . 
The unit of x is then equal to 3.6 yum, and a unit of ^ to 1.4 ms. 



case shown in Fig. [Hb). The density distribution of figure 
Ob) is essentially a continuous collision of two solitons with 
a 71 phase-diflTerence, and has the appearance of two individual 
wavepackets that never cross. These are expected to be better 
described by alternative treatments, described by Gordon 1I22I] 
in the absence of any external potential, and by Gerdjikov et 
al [ 21] in the case of harmonic confinement, or may alterna- 
tively be treated by perturbative approaches ll24l[25l l26h. 

In 3D simulations, collisions of solitons with a tt phase dif- 
ference have been predicted to have some degree of immunity 
against collapse, as opposed to in-phase collisions 1 15, 40]. 
It is clear from the above analysis that systems of two soli- 
tons of equal norm will be able to maintain a tt phase diflfer- 
ence for all repeated collisions for any choice of initial posi- 
tions and momenta, provided the right initial phase difference 
is chosen, and thus remain stable. Systems of two solitons 
with diflTerent norms have not been studied here, but are an 
area for future study; we recall that the particle model is ex- 
pected to remain valid for regimes of different norms as long 
as 21^71 - 7/2 1 <^ |vi - V2I, as explained in Sec. lIIIB 2l 



C. Three solitons 

1. Particle model 

Whereas for two solitons, the particle model dynamics are 
always regular (see Sec. IIIIB 21) , in the case of three solitons 
(A^^ = 3), the situation is quite different. A useful coordinate 
system for the three soliton system is to be found in the nor- 
mal coordinates of the system for small displacements of the 
particles from the origin: the centre-of-mass position 

^ ^1^1 + ^2^2 + ^3^3 

Z.T := , (l^j 

771 + 7/2 + 7/3 

7/1(7/2 + 27/3)^1 + 7/2(7/3 - 7/0^2 - ^3(^2 + 27/0^3 

Zc •= — , (20) 

(corresponding to the "stretch" mode), and 

Zr := qi - 2q2 + ^3, (21) 
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FIG. 7: (colour online). Trajectories in the particle model (lines) plotted over density distributions predicted by ID GPE dynamics, corre- 
sponding to the trajectory marked on figure Hd). The relative phase of the solitons in the wave dynamics is zero in figure (a), and n in figure 
(b). The solitons have equal eff'ective masses, the axial trapping frequency is 10/2;r Hz, and the other parameters (radial trap frequency of 
800/2;7r Hz, atomic species mass and scattering length of ^Li, and 5000 particles per soliton) are comparable to those in recent experiment [ 13ll . 
The unit of x is then equal to 3.6 jjm, and a unit of ^ to 1.4 ms. 

(corresponding to the "asymmetric stretch"). The stretch to ID, however, there is no analogue of the molecular bending 
modes are similar to those used to describe vibrational dynam- mode. Using these coordinates, the three-particle Hamiltonian 
ics in a tri-atomic molecule f66^ ; as the system is constrained [Eq. ([T6b l takes the form: 



"=2 



■ + wt 



■ + w 



2mm + mm + ^mm 



OJ 



m+m + m "^mm + mm + ^mm ' mmm 

'^■^^2, . .X. imm-^mm-^^mm . 2 mmm 

ZT(m-^m-^m)-^Zc — + — — 

7/1 + 7/2 + 7]3 r]i7]2 + mm + 4^1^3 _ 

^mm I mm + 2^71^3 



- 2r]iT]2(r]i + 7/2)sech^ 

- 2t]iT]3(t]i -h 7/3)sech^ 

- 2r]2r]3(m + ^3)sech^ 



m +m\mm + mm + 4^71^2 
2riim I mm - ^mm 



m +m\mm + mm + ^mm 
2mm I -(mm + 2mm) 



m + m\mm + mm + 4^71^2 
I 



Zr+Zc 



Zr + 2Zc 



Zr+Zc 



(22) 



where Wt = pi+ P2 + P?>,'^c = {(Jli + 2m)Pi + im ~ m)P2 - ^limP?) I imm + ^2^3 + 4771773) are the momenta canonically 

(772 + 2m)P3]/(m + ^2 + ^3), and = (mmpi - 2mmP2 + 
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the centre-of-mass degrees of freedom removed, becomes: 
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FIG. 8: (colour online) (a) Difference between soliton peak trajectory 
from wave simulation and particle model. The trajectories are those 
of the left-hand soliton in figure |2] going into and re-emerging from 
a collision in the case of A(p = tt (dark line) and A0 = (light line). 
Zero is indicated by the dotted line. The diff'erence Aq is equal to the 
particle trajectory qi minus the position of the peak of the left-hand 
soliton, or equivalently the position of the peak of the right-hand 
soliton minus the particle trajectory q2. (b) Diff'erence between the 
curves in figure (a). Note that the separation between the trajectories 
in the particle model and GPE dynamics is always larger in the tt- 
phase case. 



conjugate to the coordinates Zt, Zc, and Zr, respectively. It 
is apparent in Eq. (l22l) that the Hamiltonian, as in the two- 
particle case, is decoupled into a centre-of-mass component, 
and a component describing the stretch modes (which are cou- 
pled to each other). 

In the case of identical effective masses, the coordinates 
simplify substantially. It turns out to be convenient to con- 
sider slightly different coordinates, however, as this produces 
a simpler final form for the Hamiltonian describing the stretch 
mode dynamics. We therefore define Qt = tjZt = rjiqi + ^2 + 
^3)/3, qc = rjzc = r](qi - q^)l2 and qr = 7/Zr = niqi + ^3 - 2^2). 
We rescale the time to f = rf't, and then introduce the mo- 
menta pc = Wclrf- = {pi - p^)lrf' and pr = Wrlrf = 
{pi - 2p2 + P3)/6t]^. Using these dynamical variables, the re- 



2 2 



pI 



- 4sech2(^, + |) - 4sech2(^, - |). 



(23) 



This Hamiltonian, describing the two remaining degrees of 
freedom, is not separable, and it is necessary to integrate the 
corresponding Hamilton's equations of motion numerically 
to analyse the system's behaviour. As they represent a slice 
through the phase space of a system, Poincare sections pro- 
vide a good illustration of regions of regular and chaotic dy- 
namics. In regions of regular behaviour, any trajectory will lie 
on a torus in phase-space, and will thus trace a closed curve in 
the Poincare section; in regions of chaotic behaviour, a trajec- 
tory will go through every point in that region of phase space, 
and thus fill an area on the Poincare section (a so-called er- 
godic sea) |64, 65]. We choose to show sections correspond- 
ing to the momentum pr and position qr of the "asymmetric 
stretch" mode when the "stretch" mode coordinate takes the 
value qc = 0, and when its canonically conjugate momentum 
Pc < 0. Other sections can be expected to be equally illustra- 
tive of the qualitative behaviour. 

Figure [9] shows six Poincare sections for three different re- 
duced system energies H. The behaviour is regular at large 
positive values of H [Fig. |9](a)], but as H is reduced, chaotic 
behaviour emerges, characterised by ergodic regions in be- 
tween regular tori. For small (negative) H the system is mostly 
an ergodic sea, with islands of stability [Fig. |9](c)]; but as H 
is made more negative, the chaotic regions begin to subside, 
and the behaviour becomes increasingly regular again. 

Consideration of the form of the reduced- system Hamilto- 
nian [Eq. (|23]) 1 shows that without the interaction the system 
is integrable, as it becomes a decoupled pair of harmonic os- 
cillators. When H is large and positive, the interaction part of 
the Hamiltonian (which is always negative) should give a rel- 
atively small contribution to the Hamiltonian, compared to the 
integrable part of the Hamiltonian (which is always positive). 
When H is reduced, this is no longer the case, and chaotic 
dynamics are manifest. However, in regimes where the co- 
ordinates and momenta are close to zero, i.e., H approaches 
its lower bound of -12, the interaction potential becomes ap- 
proximately harmonic. The Hamiltonian H takes the follow- 
ing seperable form: 



H 



--3pj + 



247/4 



+ 2 



El 

4 



27/4 



-h24 



(24) 



sultant Hamiltonian (the reduced system Hamiltonian), with 



i.e., it again describes a pair of decoupled harmonic oscilla- 
tors. We consequently expect the phase- space structure to be 
qualitatively similar in the opposing limits of H very large and 
positive, and H large and negative. From Figs. [3a) andl^f). 
we do indeed observe this to be the case. 



2. Comparison with GPE simulations 

Figure [TOl shows a comparison of trajectories in the particle 
model with results from integrations of the ID GPE [Eq. ^] 
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FIG. 9: (Colour online). Poincare section of the three-soliton system with (a) H - 60; (b) H - 10, regions corresponding to trajectories in 
figures 2(a) to 2(c) are labeled, and highlighted using larger, darker points; (c) Poincare section of the system with H - 2\ {&) H - -2; (e) 
H = -5\ {f) H - -\0 The section corresponds to the momentum pr and position qr of the "asymmetric stretch" mode when the "stretch" mode 
coordinates qc -Q,Pc < 0. The figures correspond to the regime where the solitons have equal eff'ective masses, the axial trapping frequency is 
10/2:71 Hz, and the other parameters (radial trap frequency of 800/2;7r Hz, atomic species mass and scattering length of ^Li, and 5000 particles 
per soliton) correspond to recent experiment 1 13]. 



for the three-soliton system, where the solitons all have equal 
effective masses. As with the two-soliton case (Sec. lIIIBl) , the 
trajectories in the particle model gradually acquire a shift with 
respect to the trajectories traced out by the GPE wavef unction 
peaks. In Figs.[TOta) andfTOtb) the overall shift indicates that 



subsequent collisions tend to take place sooner in the particle 
model than is predicted by the GPE evolution; interestingly, 
in Fig. [Tote) we observe the opposite, however. As before, 
these shifts are caused by an accumulation of small errors, 
due to the fact that within a harmonic confining potential the 
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FIG. 10: (colour online). Trajectories in the particle model (lines) plotted over density distributions predicted by ID GPE dynamics, corre- 
sponding to a regular orbit. The parameters of the system are ^=10, the solitons have equal effective masses, the axial trapping frequency is 
\0/27T Hz, and other parameters (radial trap frequency of 800/2;7r Hz, atomic mass and scattering length of ^Li, and 5000 particles per soliton) 
correspond to the recent experiment (.13.1 . The unit of x is then equal to 2.4 //m, and a unit of t to 0.6 ms. 



individual solitons do not move asymptotically far from each 
other subsequent to collisions. In Fig.[TOtc) there is the added 
complication that two of the solitons appear to have formed a 
"bound state." 

The comparisons illustrate the good agreement between the 
particle model and the ID GPE in the regimes in which the 
particle model is valid, i.e., when solitons are well separated 
between collisions [Fig. [TOta) andfTQtb)1, even when the mo- 
tion is chaotic [Fig. [TQtb)1. When two of the solitons are not 
well separated [Fig. [TOtc)], the ID GPE simulation shows that 
a "bound state" is formed, which looks like a single "higher- 
order" soliton with an excited breathing mode |23]. The par- 
ticle model does not predict well the behaviour within the 
"bound state", but does give a good prediction of the centre- 
of-mass motion of the "bound state" and its interactions with 
the other soliton; it is likely that the behaviour of the den- 
sity of the "bound state" is strongly coupled to the phase be- 
haviour within the "bound state." As in the two soliton case 
(Sec, nil Bt , errors gradually accumulate in the particle model 
which lead to an overall "time shift" in the overall collision 
dynamics. It should be noted, however, that apart from this 
shift qualitative agreement with the dynamics predicted by the 
GPE remains quite good right up until limits of our numerical 



calculations (t = 10000 corresponding to 6 seconds for the 
experimental parameters in |[T3ll ). 



IV. CONCLUSIONS 

In this paper we have investigated the soliton-like nature of 
two and three bright solitary waves in harmonically trapped 
ID Bose-Einstein condensates. To this end we employed a 
model treating each soliton as a classical particle and com- 
pared the particle trajectories with simulations of the Gross- 
Pitaevskii equation. The results from the particle model and 
the GPE display good agreement when the solitons collide 
with large relative velocities, such that the collision time is 
small with respect to the trap period. When the solitons' 
relative velocities are reduced, the trajectories in the particle 
model "get ahead" of those in the GPE simulations, i.e., points 
on the particle trajectories can be identified with correspond- 
ing points in time during the GPE evolutions, but the GPE 
evolution is increasingly delayed by comparison; this eff'ect is 
more pronounced when the solitons collide out of phase, and 
is due to the non-zero time in which the particle trajectories 
asymptote to the wave trajectories after collisions. 
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We have shown that, when the external potential is har- 
monic, for systems of two solitons of equal size, the phase 
difference between the solitons is preserved for repeated col- 
lisions. In these systems, repeated out-of-phase collisions can 
cause the discrepancy between the dynamics in the particle 
model and the GPE to build up relatively rapidly. We note that 
it is in this regime, i.e., solitons colliding n out of phase, that 
the solitons are predicted to be stable against collapse when 
3D effects are considered. 

Finally, we have extended the treatment to systems of three 
harmonically trapped solitons. Using the particle model we 
have shown that, unlike the two soliton systems, systems of 
three solitons can display chaotic dynamics. Chaotic regimes 
rescind when the energy of the particle system (with the 
centre-of-mass motion neglected) is very positive or very neg- 
ative. In these limits, the dynamics decouple to those of two 
separated simple harmonic oscillators. The agreement be- 
tween the dynamics in the particle model and the GPE sim- 
ulations is good in both regular and chaotic regimes. 



where 



%j — ( 



^iD = g3D dydz\^{ymz)V = 2ncDra. (A5) 



Reassigning the zero of energy allows us to drop the con- 
stant term ficoy from the Hamiltonian [Eq. (IA4I) 1 . Introducing 
the dimensionless variables 



and 



m\giY)\N 



ft3 



(A6) 



(A7) 



Eq. (IA4b then becomes: 

= ^ + Y^¥« - l^«l¥«. (A8) 
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APPENDIX A: FROM 3D TO ID GROSS-PITAEVSKII 
EQUATION 

The one-dimensional GPE can be obtained in approxima- 
tion from the three-dimensional GPE by assuming a gaus- 
sian ansatz for the radial wavefunction and integrating out the 
radial dimensions. The gaussian ansatz can be taken to be 
the approximate radial solution if the radial potential is suf- 
ficiently tight that the harmonic potential energy dominates 
over the interaction energy in the radial directions ll4lll42ll . 

Beginning with the 3D GPE: 



ot 



2m 



V' + yext(r) + g3DWr)P 



where 



*P(r) (Al) 



(A2) 



and *P(r) must have unit norm. We employ the ansatz *P(r) = 
ip{x)(S>{y)(S>{z), where 



1/4 



exp 



20-2 



(A3) 



is the harmonic oscillator ground state, with cr^ = 
h/mojr. Averaging over the radial degrees of freedom 
£^ dydz<^*(y)<^*(z), Eq.lATIbecomes: 



ot 



2m dx^ 



2 2 

mojzx 



+ giDN\i//(x)\^ + najr 



(A4) 



m := —=^(x), 
^JmgN 



(AlO) 



such that iff(x) is normalised to one with respect to x. 

In the case where o) = 0, Eq. (IA8I) has the general soliton 
solution which takes the form, in the limit t ^ ±oo |JJ, 

H^) = ^{27/ysech [27]j(x -vjt- xoj)] 

j 

X exp [ivj(x - Vjt- xoj) -h (2r]j -h vj/2)t-\- iaoj]} 

(All) 

The normalisation requirement assumed in Eq. dASb implies 
that 2; ^7 = 1/4. 

Throughout the body of the paper, whilst the ID system is 
being discussed, the tildes have been dropped for notational 
convenience. 



APPENDIX B: HARMONICALLY OSCILLATING 
SOLUTION OF THE ID GPE 

Here we show that an arbitrary solution to the ID GPE with 
harmonic potential 



i//(x, t) = u(x, t) exp [/0(x, 0] 



(Bl) 



can be converted to an oscillating solution with a modified 
phase: 

</^(x, = u{x - (x), t) exp {/ [0(x, -h p(t)x - S (t)]} . (B2) 
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This general result is used twice in this paper; firstly to infer 
the phase behaviour of a single soliton. Secondly it is used 
to show that any trajectory of two identical solitons is identi- 
cal, in its collisional density, to the corresponding symmetrical 
case in which the solitons repeatedly collide in the centre of 
the potential. We note in passing that this result holds for any 
nonlinear Schrodinger equation with the nonlinearity given by 
some function of only. 

If ij/^x, t) = u(x, t) exp 0] is a particular solution to the 
ID GPE, where u{x, t) and t) are real- valued functions, 
the following coupled equations describe their behaviour: 

du 



dt 



dx dx dx^ 



'di 



dx 



-h -H(x)u 
u 



(B3) 



(B4) 



Defining the functions u(x,t) := u(x -h <x), and 6(xj) := 
(p(x-\-{x}), where (x) = xq cos(ojt) -\- (po/oj) sin(6i;0, and defin- 
ing the coordinate ^ = x - (x) the following relations follow 
trivially: 



u(^, t) = u(x, t), 
6(^j) = cp(x,t), 
du(x, t) du(^, t) 



dx 



du(x, t) ■ ^ du(^, t) du(^, t) 
= -{x} r— + 



dt 



dt 



(B5) 
(B6) 

(B7) 
(B8) 



Also, the relationship between the partial derivatives of ^(^, t) 
and t) is the same as that between t) and u{x, t). From 
now on it is assumed that u, and their derivatives are be 
functions of ^ and t. 

Eqs. (|B3]) and (lB4l) become: 



du 
'dt 



du I do 
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do 1 Ide 



dt \d^ 
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S<^>-¥-lSj -.Hm-^ixf-co^m. (BIO) 

By choosing ^(^, t) such that 



and 



where 
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S{t) = \pl- 
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sin(2a;0 x^p^ xopo 

— -\- — — cos(2a;0 - — — -\- S (0), 

zoj 2 2 

(B13) 

we can relabel ^ as x, and find that (p(x, t) and u{x, t) are solu- 
tions of Eqs. (IB3I) and (IB4I) . Hence, u{x, t) has the profile of 
u{x, t), but undergoes additional global harmonic oscillations 
at the trap frequency, and the result is proved. 



APPENDIX C: MOTION OF A SINGLE TRAPPED SOLITON 



The motion of a single soliton can be derived by adapting 
the method of Scharf and Bishop |32] where a solution to the 
homogeneous GPE is used as an ansatz for the GPE with an 
external potential. This method has the advantage that it can 
be used with other than harmonic potentials. Here we show 
that for a harmonic potential it confirms the expected result of 
simple harmonic motion. 

In the ID GPE, the norm TV, and energy 6, are conserved 
quantities. The single soliton solution of the homogeneous 
case [Eq. (|7]) with j = 1] is used as an ansatz for the harmonic 
case. Evaluating the norm and energy functionals with this 
ansatz yields: 
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(C2) 



The conservation of the norm N, leads to 7] ^constant, and the 
conservation of the energy 6, leads to an equation of motion 
for the peak of the soliton, q: 



q = -oj q. 



(C3) 



APPENDIX D: DEDUCTION OF THE INTER-PARTICLE 
POTENTIAL 

In order to deduce an eff'ective inter-particle potential, it 
is necessary to perform a classical inverse- scattering calcula- 
tion to find the potential which produces the same asymptotic 
position-shifts as given in Ref. 1 5] for solitons emerging from 
a collision. This proof follows the method given in Ref. [621 . 

The position-shifts given in Ref. |5] are equivalent to the 
asymptotic time- shifts for two solitons of initial speeds vi and 
V2, and eff'ective masses 7]i and 7/2, given by 



1 



2(vi - V2) 



(1.1) 

ml 



In 



(vi 



■V2f +4i7]i+7]2f 



(Vl 



V2f + 4 (r/i - 7]2f 



(Dl) 

We wish to produce these shifts in a system of two classical 
particles described by the Hamiltonian: 



H 



2fir 



+ V{q) 



(D2) 



where q := qi - q2 is the relative coordinate,//^ := ^1^2/(^1 + 
7/2) is the reduced eff'ective mass, p = (7]ip2 - ^2/^1)/ (^1 + ^2) 
is the relative momentum, and the centre of mass has been 
separated from the problem. For particles initially separated 
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at infinity, and noting that = this Hamiltonian 

takes the asymptotic form 



(D3) 



i.e., we assume the potential must vanish asymptotically. 
By rearranging Eq. (ID2I) . we may write the infinitesimal 



dq 



\2[Eoo-V{q)] 



(D4) 



since energy is conserved over the whole trajectory. Integrat- 
ing the time diff'erence between trajectories with and with- 
out the inter-particle potential, we determine the asymptotic 
timeshift to be 



-1/2 



- 1 



(D5) 



Now, expanding Eqs. (IDlb and (ID5t in terms of powers of 
l/Eco, and equating equal powers, we obtain 



r dqViqT =U- ^ -] 
J-oo 2\r]i 772/ 



X I [47/1 7/2 (^1+^2)]"- 

(-ir2«-i [(n-l)lf 
^ (2^-1)! 



47/17/2 (7/1 - 7/2) 



(D6) 



for all positive integers n. We now evaluate the integral of a 
candidate potential: 



2 / 27/17/2 

2(7/1 +7/2)7/1 7/2 sech ' 



dqV{qr= dq\- 

=^(- + -)[4^i^2(^i+^2)r 

2\^i ^2/ 



(272-1)! 



(D7) 



Comparing the expressions (ID6b and (ID7I) it follows that the 



potential 

V{q) = -2(7/1 + 7/2)7/i7/2sech2 (^^^) (D8) 
gives the correct time shift in the limit 2|7/i - 7/2I «: |vi - V2I. 



APPENDIX E: COLLISIONAL FORM PRESERVATION 
WITH TWO HARMONICALLY TRAPPED SOLITONS 



Let us consider a general parity operator P(0), such that 
P((p)i//(x, t) = e^^il/{-x, t) := x(x, t). We want to know whether, 
if ;^(x, ^o) = <A(-^, ^0), it continues to be the case that;^(x, t) = 
i//(x, t) for all t. 

We take the time derivative of ;^(x, 0- Noting 
that P{(p)\il/{xj)\^il/{xj) = e''^\i//(-x,t)\^i//(-xj) 
\P((p)i/^(xj)\^[P((p)i/^(xj)] we deduce from the ID Gross- 
Pitaevskii equation [Eq. ©] that 



d d 
i^xixj) =iP((f>)—i//(x,t) 
at at 

1 52 



P{<P) 
"2 5x2 



2 2 

CO X 



2 5x2 2 



2 2 



X(x, t). 



(El) 



Hence, we see that the time-evolutions of i/^(x,t) and 
P((p)i//(x, t) = x(x^ are governed by the same diff'erential 
equation. If we also choose an initial condition such that 
i//(x, to) = x(x, to), it must therefore follow that;^(x, t) = i//(x, t) 
for all t. In other words, parity is conserved in the sense that 
an initially symmetric wave function will have that symmetry 
preserved throughout its subsequent dynamical evolution. 

In this paper, the most important consequence of this result 
is that a system of two identical solitons with equal and op- 
posite velocities will repeatedly collide with the exact same 
collisional form (e.g., in phase or tt out of phase) at the exact 
centre of the trapping potential. Using the results of appendix 
IHl it follows that an equivalent result holds upon the addition 
of a centre of mass oscillation. 
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